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The universal traveling wave solution to the Balitsky-Kovchegov equation with running coupling 
(and other equations in the same universality class) is extended to subleading orders at large rapidity 
and small dipole size r. The large rapidity expansion of the logarithm of the saturation scale Q a (Y) 
is derived from that traveling wave solution. In addition to the two already known leading terms 
in y 1/2 and Y 1/6 , which are determined only by the LL BFKL kernel, the three following terms in 
Y°, Y~ 1/6 and Y~ 1/:i are obtained. They are universal and sensitive to NLL BFKL effects. Initial 
condition dependence and NNLL BFKL effects would both start to appear only at the next order, 
which is in Y~ 1//2 . In the light of these results, the impact of the precise implementation of the 
running coupling in gluon saturation is discussed, i.e. parent dipole size prescription vs. Balitsky's 
prescription, and also the effect of the full NLL corrections, not yet implemented in numerical 
simulations. 
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CLf I. INTRODUCTION 

or 

In the high energy limit for hadronic collisions, softer and softer partons are resolved in the hadronic wave- 
. functions. According to the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [1-3], these soft partons are emitted 
in a branching process [4]. Eventually, the occupation number of the soft gluons becomes large, leading to the 
phenomenon of gluon saturation [5-9] (for recent reviews see e.g. [10-13]). Further emission of softer gluons is a 
J> ' nonlinear collective effect, described by the Balitsky-Kovchegov (BK) [14-16] or JIMWLK 1 [17-24] equations, which 
generalize the BFKL evolution and overcome two of its drawbacks, i.e. the violation of unitarity at fixed impact 
parameter and the diffusion into the infrared. The boundary between the saturated part of a hadron wave-function, 
consisting of partons with small transverse momentum k± and/or large rapidity separation Y and the dilute part is 
parameterized by the saturation scale Q S (Y). 

Some properties of the solutions of the gluon saturation equations such as BK and JIMWLK can be derived 
analytically [5, 25-28]. The most powerful method for this purpose is borrowed [29] from the studies of nonlinear 
wave-front formation in some reaction-diffusion systems [30, 31]. The most interesting result is that many features 
of the solutions of the saturation equations are universal, i. e. independent of the initial condition, such as the large 
Y behavior of the saturation scale, in agreement with numerical simulations [32-34] . 

Early phenomenological studies have found gluon saturation qualitatively consistent with the deep inelastic 
scattering (DIS) data at low XBj from HERA. For example, gluon saturation naturally explains the geometric scaling 
property of the data [35]. Later, more obvious indications of gluon saturation have been found in particle production 
at forward rapidity in d-Au collisions at RHIC, such as the suppression of the Cronin peak [36, 37] or recently the 
disappearance of the recoil jet [38] in two particles azimuthal correlations. However, the phenomenological rapidity 
dependance of Q S (Y) extracted from the data is much slower than the one predicted from the BK and JIMWLK 
equations at leading logarithmic (LL) order. 

That observation suggests to consider the effect higher order corrections to these saturation equations. Next-to- 
leading logarithmic (NLL) corrections to the BFKL equation [39, 40] are known and indeed large. The quark part 
of NLL corrections to the BK and JIMWLK equations has been obtained a few years ago [41-43]. Most of these 
corrections correspond to running coupling effects. One can then build improved LL BK and JIMWLK equations 
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with running coupling at an appropriate scale [41, 42]. Recently, the full NLL corrections to the BK equation have 
been derived [44, 45] . State of the art phenomenological studies are now based on numerical simulations [46] of the 
improved LL BK equation with Balitsky's running coupling prescription [42], and provide a unified description of the 
DIS data [47] , single inclusive particle production [48] and two-particles azimuthal correlations [49] in d- Au collisions 
at RHIC, and two-particles long range rapidity correlations in Au-Au collisions at RHIC [50]. 

Analytical results are important to gain some insight into these numerical studies, in particular to understand 
the potential impact of the full NLL effects besides the running of the coupling, which are not yet implemented 
numerically. Previous studies [5, 25-27] have shown that in the running coupling case, the solutions of the saturation 
equations seem still universal, but are in a different universality class than the fixed coupling ones, and the two leading 
terms in \ogQ s (Y) at large rapidity have been obtained. It has been shown [51, 52] that these two leading terms are 
independent of the details of the running coupling prescription and of the NLL (or higher order) contributions not 
associated with the running of the coupling. Hence, one has to calculate further sublcading terms in the asymptotic 
behavior of the solutions and of logQ s (y) in order to discuss full NLL effects in a consistent way. That observation 
raises some questions about the validity of the previous attempt at obtaining NLL effects on the saturation scale 
evolution [53]. 

In this letter, the results for three new universal terms (19,20,21) in the asymptotic expansion (18) of \ogQ s (Y) in 
the running coupling case are presented in section IV and their derivation is outlined. These new terms are NNLL 
independent, but two of them are NLL dependant. The corresponding asymptotic expression for the saturation scale 
is used in section V to compare the outcome of various saturation equations: the LL BK equation with the parent 
dipole size prescription or with Balitsky's prescription for the running coupling, and the NLL BFKL equation in 
momentum space supplemented by saturation effects. But first, let us start in the next section by a review of the 
known results in the fixed coupling case, as a warm-up. 



II. REVIEW OF THE FIXED COUPLING CASE 



At leading logarithmic (LL) order, the gluon saturation equations such as BK or JIMWLK can be written formally 
as 2 

dyN(L,Y) = a { x (-9l)N(L, Y) - Nonlinear damping} (1) 

with the notation a = a s N c /n 7 both in position space or momentum space. In the former case, N(L, Y) is the dipole- 
target amplitude for a dipole of size r and a rapidity interval Y, with L = — log(r 2 Q 2 /4) and Qa an arbitrary reference 
scale. In the latter case, N(L, Y) is typically an unintegrated distribution of gluons with transverse momentum k± 
and rapidity Y in the target, with L = log(k± 2 /Qq 2 ). The linear part of equation (1) involves the BFKL kernel at 
LL accuracy and restricted to zero conformal spin, whose characteristic function is 

X (7) = 2*(l)-*(7)-*(l-7) (2) 
in both the position and momentum spaces, ^(7) being the digamma function. 

The BFKL linear term in (1) implies an exponential instability in Y of the vacuum solution N(L, Y) = and 
a diffusive behavior in L-space. The non-linear damping term in equation (1) tames the instability when N{L,Y) 
becomes of order 1. More precisely, we have N(L,Y) < 1 for the dipole-target amplitude in position space. For 
its Fourier transform as defined in Ref. [16], the exponential growth becomes a linear growth in Y in the deep 
saturation regime. By contrast, the effective unintegrated gluon distribution appearing in the fc^-factorization for 
gluon production in DIS and pA turns over in the nonlinear regime and decreases back towards zero [59]. For 
any of these cases, one expects a perturbative power tail in kj_ or r in the UV, meaning a largc-L exponential 
tail of N(L,Y). Hence the nonlinear regime of saturation is reached in the IR at lower values of L first. At 



2 For simplicity, let us assume here that the target is homogeneous is the transverse plane. However, the traveling-wave method discussed 
here is valid in the general case of scattering at non-zero transfer on an inhomogencous target[54, 55]. In this paper, the gluon saturation 
equations including gluon number fluctuations are not considered. These fluctuations occurring in the dilute regime have a big impact 
at fixed coupling [56]. When the coupling is running, these fluctuations should have similar effects in the high rapidity limit [57, 58], 
but arc completely negligible in the relevant rapidity range. 
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large enough finite values of Y, the nonlinear term in equation (f) is thus negligible for large L but not for 
small L. The typical boundary value between these two regimes L S {Y) = \og(Q s (Y) /Qq 2 ) defines the satura- 
tion scale Q S (Y). Due to the exponential rise of N(L,Y) with Y at large L, Q S (Y) is expected to grow with 
the rapidity Y. Hence, the relevant solution for N(L,Y) is a traveling wave front in L-space, where aY plays 
the role of the time variable. At the time aY, the wave front is at the position L = L S {Y) and moves towards larger L. 

As noticed in [29], the properties of exponential instability and diffusion in the linear regime and of non-linear 
damping are shared by both the equation (1) and the Fisher and Kolmogorov-Petrovky-Piskounov (FKPP) equation 
[60, 61], which can actually be written as (1) with x(l) = 1 + 7 2 and a quadratic or a cubic nonlinear damping. Hence 
some mathematical results [30, 31] about nonlinear wave front formation for the FKPP equation remain valid in general 
for equation (1). Both of these equations admit a special family of exact scaling solutions N(L,Y) — N 7 (s v (L,Y)), 
also called uniformly translating front solutions, with 

the scaling variable s v (L, Y) = L — vaY , (3) 
a scaling function with a tail N^(s v ) cx e~ 7S " for large positive s v , (4) 

y('y) 

and the dispersion relation ^(7) = . (5) 

7 

An important feature is that the wave- front velocity ^(7) admits a positive global minimum v c = v(-y c ) in the relevant 
range for 7, which is < 7 in the FKPP case and < 7 < 1 when using the BFKL eigenvalue (2). The scaling solution 
N 7c (s Vc (L,Y)) with minimal velocity is called the critical solution. The minimization condition for the velocity 

Xilc) =lcX\lc) (6) 
gives 7c = 1 and v c = 2 in the FKPP case and -f c ~ 0.6275 and v c ~ 4.883 in the QCD case. 

The generic solutions of the FKPP and BK equations, outside that family of uniformly translating fronts, fall into 
two classes [30, 31]. 

• If the initial condition has a flatter tail than the critical solution at large L, then the long Y behavior of N(L, Y) 
is very sensitive to the initial condition. 

• If the initial condition has a steeper tail than the critical solution, N(L, Y) converge to the critical front solution 
at large Y, and is called a pulled front solution. Moreover, the convergence happens first in the fully nonlinear 
regime N(L,Y) = 0(1), and propagates in an universal way into the tail via diffusion. Hence, the dynamics 
of the solution is determined by the existence of the nonlinear term in equation (1) even in a wide kinematical 
region where its magnitude is actually negligible compared to the magnitude of the linear terms. 

In QCD, the tail of the initial condition should be of the type exp(— L), due to the property of color transparency, which 
is steeper than the critical tail exp(— ~/ c L). The physically relevant solutions for QCD are thus of the pulled front type. 

In Rcf. [31], a method has been developed in order to study universal properties of the pulled front solutions of 
equations in the universality class of the FKPP equation and of equation (1). Large Y asymptotic expansions of the 
pulled front solutions can be constructed in two different regimes. The first one is the front interior regime where 
L—L S (Y) is kept fixed when Y — > 00, corresponding to describe the wave- front in its rest frame. The second one is the 
leading edge regime where L~L S (Y) = 0(V aY) when Y — > 00, corresponding to the crossover between the universal 
region of the front and the part of the tail still driven by the initial condition. Matching the two expansions together 
and imposing some causality requirement in the leading edge regime allow to determine most of the parameters and 
extract universal information [31]. For example, N(L,Y) depends only on the scaling variable L — L S (Y) in the 
increasing window L — L S (Y) <C y/2 x"(lc) aY . Combined with the appropriate dipole or fc_L-factorization formula, 
gluon saturation thus explains the geometric scaling property seen in the low- a; HERA data [35] . One also obtain the 
universal asymptotic expansion of the position of the front [31], related to the saturation scale, 

Within QCD context, the first term in the expansion (7) was derived in [5, 25], the second one in [26, 27] and last 
universal term in [28]. Due to the L-translational invariance of equation (1), traducing the conformal invariance of 
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high energy QCD at LL accuracy, the constant term in (7) is directly sensitive to the initial condition, instead of 
being universal. Moreover, equation (1) is also invariant under ^-translations, and such translations modify the large 
Y expansion of L S (Y), starting at the order 0(aY)^ 1 . Hence, in the expansion (7), the coefficients of the terms of 
order or further cannot be universal. 

In Ref.[62], the effect of NLL corrections to the kernel of equation (1) has been studied, keeping the coupling fixed. 
It amounts to deform the dispersion relation (5), and thus to shift the critical parameters -f c and v c by corrections 
proportional to a. Hence, all the coefficients in the expansion (7) are shifted each time one goes to the next logarithmic 
order in perturbation theory, LL, NLL, NNLL, etc. That property is true for gluon saturation in N=4 SYM theory. 
However, as we will see in the rest of this paper, this is not the case in full QCD, where running coupling effects 
stabilize the asymptotic expansion of log Q s with respect to other higher order corrections. 



III. SATURATION EQUATIONS WITH RUNNING COUPLING 

The most important correction in QCD to the gluon saturation equations of the type (1) is the running of the 
coupling a. That effect corresponds an all-order resummation of a tower of terms in the perturbative expansion. 
However, as a first step, it is customary to simply promote the coupling a in equation (I) to the one-loop expression 

I UN c -2N f 
a^-, where b = — , (8) 

without trying to calculate this tower of perturbative corrections. In momentum space, that prescription corresponds 
to take the running coupling at the scale of the parent gluon k±. By contrast, when applying the prescription (8) to 
an equation (1) in position space, such as the BK equation, one obtains the running coupling at a scale related to the 
parent dipole size r. By consistency, in the running coupling case one should take the reference scale Qo to be Aqcd 
in the appropriate renormalization scheme, so that now 



k± \ T fr 2 A, 



L = lo S (- r ^) or L = -!o,(-^^) . ,0) 



QCD 

The gluon saturation equations with this type of running coupling prescription have been intensively studied 
numerically, and are suitable for our analytical calculations. 

It is possible to write in a formal way the all order generalization of the gluon saturation equations (1), following 
this prescription, as an asymptotic series in (bL)~ n 

d Y N(L,Y) = ^- {x(-d L )N(L,Y) - Nonlinear damping | 

+ {^"(^W'^) _ Nonlinear dampingj 

+ .... (10) 



The usual form of the improved LL BK equation with Balitsky's running coupling prescription [42] is difficult to 
handle directly with our analytical methods. However, it is possible to expand the running couplings in order to 
rewrite that equation as a series of the type (10) [52], but this time with the NLL contribution 



Xnll{i) = ^ 



xii) 2 -x'{i)--x{i) 



(11) 



Since only part of the NLL corrections are included in that case, the renormalization scheme is not really fixed. 
However, it is natural to rescale Aqcd with respect its MS scheme value as 



Aqcd = e y ' A QCD . 
This choice corresponds to the position space analog of the usual MS scheme. 



(12) 
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Due to some specific technical issues, the full NLL BK equation [44] will not be considered here, and its study 
is postponed to further works. Let us consider instead the NLL BFKL equation [39] in momentum space for the 
unintegrated gluon distribution, and add gluon saturation effects to obtain an equation of the type (10). In that case, 
the NLL eigenvalue xnll(i) writes 3 



XNLdl) = ix(7)x'(7)-^x"(7) + ^C(3)-0(7)-^[x(7) 2 -x'(7)] + 



56+1 



7T 

12 



X(7) 



7T 2 C0S(7T7) 

4(1-27) sin 2 (7T7) 



N f \ (2 + 37(1-7)) 



N*J (3- 2 7 ) (l+2 7 ) 



4sin(7T7) 



(13) 



in the MS scheme and defining the rapidity as Y — log(s//cj_ 2 ), with the notation 



#y) = £(-ir 



n=0 



*(n+l+7)-*(l) *(n+2-7)-*(l) 



(n+7) 2 



+ 



(n + l-7) 2 



(14) 



IV. THE SATURATION SCALE FOR TRAVELING WAVES WITH RUNNING COUPLING 



The mechanism of universal pulled front formation from steep enough initial conditions has been established 
rigorously [30] in the case of the FKPP equation, and remains valid for equations within the same universality class, 
such as the BK equation at LL accuracy. By contrast, the extension to the running coupling case is far from obvious. 
There is an evidence from numerical simulations [32, 34] that some analogous universal wave front formation holds, 
but no mathematical proof. 

Qualitatively, one can justify such universality. Most of the ingredients allowing for pulled front solutions of the 
FKPP equation are still present in equation (10), such as the instability of the vacuum, the diffusion effects, the 
nonlinear damping, and also the steepness of the initial conditions relevant for QCD. A priori, the only missing 
ingredient is the family of uniformly translated front solutions N(L,Y) — N 7 (s v (L,Y)) with (3,4,5). Indeed, the 
saturation equation with running coupling (10) does not admit any exact scaling solution. In the literature [26, 27], a 
family of approximate scaling variables of the type s v (L, Y) = L—f v (Y) is usually considered, and then the methods 
developed in the fixed coupling case are applied in the running coupling case from this starting point. However, such 
choice of approximate scaling variable is equivalent, in some sense, to impose by hand the geometric scaling of the 
solution, and moreover it is not unique [63]. It is thus essential to consider also other types of approximate scaling 
variables for equation (10). On the one hand, any dependance of the final results, such as the expression for Q S (Y), 
on the choice of approximate scaling variable would refute the universality property in the running coupling case. 
Indeed, the physical Q S (Y) should not depend on arbitrary choices made when writing the asymptotic expansions 
of the solution N(L,Y). On the other hand, the geometric scaling property is often argued to be due to gluon 
saturation, and indeed arises dynamically from gluon saturation in the fixed coupling case, as discussed in section 
II. It is thus important to understand if and how such a scaling property arises dynamically in the running coupling 
case at least approximately, without imposing it by hand. 

Generically, an approximate scaling variable for the equation (10) is given by any function s v (L,Y) whose partial 
derivatives write [63] 

^ = 1 + ... and ^ = -a(L)t,[l + ...]=-^[l + ...]. (15) 

where the dots stand for terms which go to zero when L, Y — > 00 with s v (L, Y)/L — > 0. For each value of v, there is 
an infinity of such approximate scaling variables s v (L, Y). The dispersion relation (5) is found to hold independently 
of the corrective terms in (15). Hence, in the running coupling case there is a critical solution, which approximately 
scales with any variable s Vc (L, Y) of the type (15) with the same critical parameters 7 C and v c as in the fixed coupling 



3 Notice that our notations differ from the ones of Ref. [39] by exchange of 7 and 1 — 7. 
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case. 



By analogy with the fixed coupling case, the solutions with steep enough initial conditions should converge to that 
critical solution. And they should approximately scale with a variable s(L,Y) given by s Vc (L, Y) supplemented by 
subleading diffusive terms encoding the relaxation towards the critical solution. At that point, it is possible to follow 
the method of Ref. [31]. One can construct the front interior and the leading edge asymptotic expansions of those 
solutions in the dilute domain. In both of these expansions, L" 1 / 3 can be taken as the small parameter, and the 
trivial exp(— 7 C s) contribution is factored out of the solution. The front interior expansion is done at fixed value of 
s{L,Y) 1 and the coefficients are polynomials in s{L, Y). The leading edge expansion, by contrast, is done at fixed 
value of z oc s(L, F)/L 1 / 3 , and the coefficients are more complicated functions of z, involving the Airy function Ai(z) 
or its derivative, times polynomials. One rejects the possibility of having the exponentially growing Airy function 
Bi(z) appearing in the leading edge expansion in order to allow for a smooth connection with the non-universal tail 
of the solution, by analogy with the FKPP case [31]. From this requirement and the matching of the two asymptotic 
expansions, one determines most of the parameters appearing in these expansions and in the approximate scaling 
variable s(L,Y). The details of this calculation as well as the asymptotic expressions obtained for the solutions 
N(L,Y) will be reported and discussed elsewhere [64]. 



The saturation scale Q s (Y) and its logarithm 



L S (Y) = log 



Q 2 S (Y) 

^QCD 



are defined by a level line 



N(L S (Y),Y) = k 



(16) 



(17) 



where k is a given small constant e.g. n = 0.1 or n = 0.01. Using the front interior expansion constructed for N(L, Y), 
one can solve the equation (17) explicitly, and obtain the universal large- F asymptotic expansion 



Ls(Y) 



—) + —{ D —) + C0 + C -V6 



D 



2v c Y 



-1/6 



+ c_ 



1/3 



D 



2v c Y 



-1/3 



+ (V- 1 / 2 ) . (18) 



The first term in (18) was derived in [5, 25] and the second one in [26, 27], whereas the three next ones constitute the 
main result of the present study: 



co 

C-l/6 
C-l/3 



= £ 



7c 



+ N 



1 D K 3 3 2 3 
72^| _ 32 + Wc~ 8 Kz + 10^ 4 



(19) 
(20) 




+K 4 



#3 1_ 

27c 3 7c 2 



— + 6K 3 

7c 



D 



27 7 3 3 7c 



2 3K 3 11K 3 2 

~3 7c 2 4 7c 4 



3i ^5 + A [! + 3 7c^3] [A^o-7ciVi] + N 2 -DN 

2 lc 



In the previous relations £i stands for the rightmost zero of the Airy function, i.e. £i 
notations have been used: 



(21) 

2.338, and the following 



g _ X' ; (7c) _x"{lc) 
2 X(7c) 2 7c w c 



2 X ( " ) (7c) , 
nlX (7c) 



and 



" ~ n\ xile) b 



for n > . 



(22) 



The parameter D measures the diffusion rate with respect to the instability growth. Each parameter K n contains 
the nth derivative of the LL eigenvalue ^(7), so that the K^s correspond to the contributions beyond the diffusive 
approximation of the LL kernel which is sometimes performed. The parameters N n quantify the contributions of the 
NLL eigenvalue xnll{i) or its derivatives to the expansion (18). 



In the expressions (19) and (21), S is a parameter depending on the height n at which the saturation scale is defined 
via (17). It can be expressed as 



S = 



(23) 
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where W_i(x) is the —1 branch of the multivalued Lambert function W(x) [65], and A is the overall normalization 
constant of N(L,Y) appearing in the front interior and the leading edge expansions. The calculations of these two 
expansions are performed using the linearized version of the equation (10), so that A is undetermined at that stage. 
The expression of N(L, Y) obtained from the front interior expansion has a maximum and decreases in the IR. This 
maximum signals the breakdown of the front interior expansion due to the onset of the fully nonlinear regime. The 
normalization constant A should then be determined by matching numerically the front interior expansion to full 
numerical solutions extending into the deep saturation regime. Hence, A is determined by the height of N(L,Y) at 
which the nonlinear damping terms in the equation (10) become relevant. The constant A is thus independent of the 
initial condition, as long as the latter is steep enough to have convergence towards the critical solution, and is also 
independent of the presence of NLL or higher order terms in the equation, apart from the running of the coupling. 
A priori, one can expect a value of order one for A. However, this value should depend on the considered quantity: 
dipole target amplitude, effective or physical unintcgratcd gluon distribution, which have a different behavior in the 
deep saturation regime. The formula (23) is valid only for j c n <C A, i.e. when the saturation scale is defined on a 
level line (17) far enough from the deep saturation regime. 

Due to many cancelations occurring either when constructing the front interior and leading edge expansions, or 
when extracting L S (Y) from them, the result (18) does not depend on the choice of approximate scaling variable 
s(L, Y) used to write the asymptotic expansions of N(L, Y). As discussed previously, this is a non-trivial consistency 
check of our method. Hence, we have a new strong hint that in the running coupling case also, the solutions with 
steep enough initial condition converge universally towards a critical one, despite the absence of a family of uniformly 
translating front solutions. 

Translations in Y would generate contributions starting at order Y~ x / 2 in the expansion (18) out of the first term. 
Due to this remark and to the Y-translational invariance of equation (10), the coefficients of the expansion (18) cease 
to be universal at order Y -1 / 2 . By contrast, translations in L are no longer a symmetry of the equation (10) due to 
the running of the coupling, i.e. due to the conformal anomaly. Hence, all the universal terms have been derived 
in the expansion (18). It has been anticipated in Ref.[66] that, in the running coupling case, the constant term c 
should be independent of the initial condition, and that initial condition effects should be delayed until the order 
y-1/2 rp ms p rec ii c tion, which was already found consistent with numerical simulations [32, 34], is thus confirmed by 
direct calculation. In practice, that property implies that the nuclear enhancement of the saturation scale Q s by a 
factor A 1 / 6 is progressively washed out by the rapidity evolution, so that at very high energy all hadrons and nuclei 
have the same saturation scale, with a normalization set by the conformal anomaly of QCD and not by the nature 
of the target. By contrast, the nuclear enhancement of Q s would survive at high energy in a conformal Yang-Mills 
theory, hidden in the non- universal constant term of the expansion (7). 

The front interior and the leading edge expansions are asymptotic expansions in powers of L -1 / 3 . It is thus 
clear that, apart from the running of the coupling, NLL BFKL effects may start appearing in the expansion 
(18) of L S (Y) at order Y°, NNLL BFKL effects at order Y" 1 / 2 , NNNLL BFKL effects at order Y" 1 , and so on 
[51, 52]. The universal terms can thus depend only on the known LL and NLL BFKL kernels. On the contrary, 
the fact that the coefficient c_i/ 6 (20) is independent of NLL contributions is an unexpected outcome of our calculation. 

The coefficient c_!/ 3 (21) depends on NLL BFKL effects only through the combinations No—j c Ni and N2—DN0. 
It is instructive to notice that, if the NLL characteristic function xnll{i) is proportional to the LL one x(l)> 
both N — -f c Ni and N 2 — DN vanish. Different choices of renormalization scheme give characteristic functions 
Xnll{i) differing by terms proportional to x(l)- Hence, the coefficient c_!/ 3 depends on the renormalization-schcme 
independent terms only of the NLL BFKL kernel with running coupling. By contrast, the constant term Co (19) 
is renormalization-scheme dependent via No, but in such a way that this compensates renormalization-scheme 
dependence of the Aqcd scale included in the definition of L S (Y) (16). Finally, one concludes that the asymptotic 
expression of the saturation scale Q S (Y) deduced from (18) is renormalization-scheme independent. 

Presumably, phenomenological studies or numerical simulations of gluon saturation at NLL accuracy would require 
a resummation of collinear higher order contributions analog to the one done in the BFKL case [67-73]. However, 
such resummations are formally NNLL effects, so that they cannot modify the coefficients (19,20,21) of the large 
rapidity expansion (18). 
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V. APPLICATIONS 



Let us now use the formal results of the previous section to discuss the qualitative effects of running coupling 
prescriptions and other NLL contributions on the asymptotic behavior of the saturation scale. For the purpose of 
our discussion, let us take the following values for the remaining parameters: k = 0.1, A = 0.8, and Nf = 5. By 

consistency, one should then take A.qq D = 0.09 GeV, in order to have a s (M z ) = 0.1184 [74]. The evolution of Q S (Y) 
is often parameterized by 

since in the fixed coupling case, A(Y) converges fast towards a constant value. Phenomenological fits of HERA data 
with models featuring a fixed A, inspired by fixed coupling gluon saturation, give typically a value of A around 0.3 or 
0.2, depending on the details of the model [75-80]. 

The results for X(Y) and Q S (Y) obtained from the expansion (18) for improved LL BK equations with running 
coupling are presented on Fig.l, both with the parent dipole size prescription, i.e. taking Nq — N\ = N 2 = 0, and 
with Balitsky's prescription, i.e. with the coefficients N n calculated using the NLL eigenvalue (11). The rescaling 
(12) of Aqcd is done here. One can see that the choice of running coupling prescription has only a negligible impact 
on X(Y). Indeed, the coefficient cq drops from X(Y) due to the derivation. The three leading terms in the asymptotic 
expansion of X(Y) are thus independent of the details of the running coupling prescription. Moreover, there are 
accidental cancelations in the coefficient c_i/ 3 further reducing the sensitivity to the prescription: from the eigenvalue 
(11), one has N ~ -2.45, N x ~ -3.03 and N 2 ^ -19.4 but only A -7 c Ai ~ -0.55 and N 2 -DN ~ -0.027. Hence, 
the main effect of the choice of prescription is a reduction of the value of Q S (Y) by a factor exp(Ao/2) ~ 0.29 (or 
0.35 when Nf = 3) for Balitsky's prescription with respect to the the parent dipole size prescription, as one can see 
from the right plot of Fig.l. That seems to be in agreement with numerical simulations (see e.g. Fig. 4 of Ref.[46]). 

The values of X(Y) obtained with both prescriptions and with no non-universal subleading contributions to the 
expansion (18) are small enough at very large rapidity, as one can see from the left plot of Fig.l, but remain too 
large in the phenomenologically interesting range 4 < Y < 10. By tuning the non-universal subleading contributions, 
for example by adding 20 to Y in the first term of the expansion (18), it is possible to stabilize X(Y) around its 
phenomenological value down to the range of interest in Y. However, such a big shift is not very natural. Moreover, 
since the asymptotic behavior of Q S (Y) is universal, reducing the evolution rate X(Y) at intermediate rapidities 
leads to a significant increase of Q S (Y) at lower rapidities. As one can see from the right plot of Fig.l, the value 
of Q S (Y) at moderate rapidities is too high 4 when the shift is performed. Hence, even by tuning initial condition 
effects, it seems difficult to get simultaneously the correct values of both Q S (Y) and X(Y) in the relevant rapidity 
range from the improved LL BK equation with running coupling, independently of the details of the running coupling 
prescription. Such difficulty has been bypassed e.g. in Ref.[47] by multiplying Aqcd by a free parameter fitted on 
the data. 



In Fig. 2 are presented the results for X(Y) and Q S (Y) for the BFKL equation in momentum space with running 
coupling at the parent gluon k±, modified by saturation effects, both with the full NLL BFKL eigenvalue (13) or 
with the LL one only. One can see on the left hand plot that the NLL corrections on X(Y) arc sizable. They 
stabilize X(Y) down to Y ~ 15, at values around 0.2 — 0.25, even without including any non-universal corrections 
in the expansion (18). In that case, there are indeed no strong cancelation occurring in the NLL contributions 
to the coefficient c_i/ 3 : from the NLL eigenvalue (13), one gets Ao ~ —7.98, N\ ~ —3.83 and A2 ~ —128.4 and 
thus No — j c Ni ~ —5.57 and N2 — DN0 ~ —65.2. Including moderate nonuniversal subleading contributions in the 
expansion (18) by adding e.g. —2.5 to Y in the first term further improve NLL results. The values of X(Y) are now 
compatible with phenomenological results in the relevant rapidity range, where, simultaneously, Q S (Y) is reduced to 
small enough values, as we can see on Fig. 2. 

Let us summarize the results of this section. Compared to the fixed coupling case, running coupling effects signif- 
icantly improve the agreement between the theoretical predictions and the phenomenological fits for the saturation 



4 Notice that our choice k = 0.1 automatically leads to higher numerical values of Q S (Y) than the usual definitions of Q a (Y) in the 
literature. However, that effect is too small to make the upper curves of the right plot of Fig.l consistent with the expectations for 
QsiY). 
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FIG. 1: Left: \(Y) (see equation (24)) as a function of Y. Right: The saturation scale Q S (Y) as a function of Y. In both plots, 
the results for the LL BK equation with running coupling are drawn with dotted lines for the parent dipole size prescription, 
i.e. with vanishing coefficients N n , and with solid lines for Balitsky's prescription, i.e. with the coefficients iV n calculated using 
the eigenvalue (11). The upper curves on the left plot (resp. lower curves on the right plot) correspond the expression (18) 
truncated after the term of order F -1 ^ 3 . By contrast, for the lower curves on the left plot (resp. upper curves on the right 
plot), a shift Y M> Y + 20 has been performed in the leading term of the expansion (18) only. 



0,6-i 1 

'> \ vf- 




FIG. 2: Left: X(Y) (see equation (24)) as a function of Y. Right: The saturation scale Q S (Y) as a function of Y. Both plots 
are obtained from the expansion (18) applied to a unitarized version of the BFKL equation in momentum space with running 
coupling at the parent gluon k±. The full NLL BFKL eigenvalue (13) is used to calculate the coefficients N„ for the solid and 
dashed curves, whereas for the dotted curves no NLL contributions are included, i.e. the coefficients N n are taken to be zero. 
The solid and dotted curves correspond the expression (18) truncated after the term of order y -1 / 3 . For the dashed curves, a 
shift Y i — ^ Y — 2.5 has been performed in the leading term of the expansion (18) only. 
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scale and its evolution. This is expected since the wave front solutions to saturation equations with fixed or running 
coupling are in different universality classes. However, when including running coupling effects only on top of LL 
equations, there remains a tension between the theory and the fits, even when strong subleading initial condition 
effects are included. By contrast, when not only the running coupling effects but also the full NLL BFKL effects are 
included, the theoretical predictions are naturally close to the phenomenological results. In that case, only minor 
initial condition effects may be required to describe the saturation scale of a proton or of a nucleus, and its rapidity 
dependance. 

VI. CONCLUSION 

The central result of this paper is the complete universal part of the large rapidity asymptotic expansion of 
logQ s (F) (18,19,20,21), valid for evolution equations with gluon saturation and running coupling. That result is 
obtained using the traveling wave method of Ref.[31]. Among the three new terms in the expansion, two of them are 
sensitive to the details of the running coupling prescription and other NLL effects, and also to some details of the 
nonlinear damping via the initial-condition-indepcndent normalization constant A of the universal traveling wave 
solutions. Initial condition and NNLL effects would appear only at the following order in the expansion of log Q S (Y), 
which is F -1 / 2 , in agreement with previous expectations [51, 52, 66]. By contrast to the fixed coupling case, the 
normalization of the saturation scale Q S (Y) at large rapidity is thus independent on the nature of the target, and is 
instead related to Aqcd, to the NLL corrections to the BFKL kernel, and to the normalization constant A which 
has to be determined numerically. 

The details of the running coupling prescription are found to have a negligible impact on the evolution of Q S (Y), but 
an important one on its normalization, in agreement with numerical simulations [46]. By contrast, the contributions 
to the NLL BFKL kernel not related to the running of the coupling have a large impact on both the normalization of 
Q s (y) and its evolution, and they allow to obtain results in agreement with phenomenological fits without the need 
of strong initial condition effects. This remark suggests that the large rapidity expansion (18) of \ogQ s (Y) maybe 
relevant down to the rapidities relevant for collider physics. Further studies are however necessary to fully understand 
its validity range. The results presented here should also provide a useful constraint on high energy extrapolations of 
models of hadronic collisions for ultra high energy cosmic rays. 
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